NoahProject
Evan Doughty
3/22/2019
Data Set-up
#############################################################################################
setwd("~/Dropbox/CalcanealGearRatio")
calcmeasure <- read.csv("~/Dropbox/CalcanealGearRatio/LACM_CalcaneusMeasurements_2018_8_29.csv", stringsAsFactors=FALSE)
calcmeasure <- calcmeasure[calcmeasure["Measurement.Type.Manual.Photo."] == "Photo",1:14]
calcSpecies <- unique(paste(calcmeasure$Genus, calcmeasure$Species,sep="_"))
MOM.Mat <- read.csv("~/Dropbox/CalcanealGearRatio/MOM_v4.1.csv")
#remove empty rows
MOM.Mat <- MOM.Mat[!MOM.Mat$Species == "",]
#aggregate to get single species entries
MOM.Mat[,"TaxonomicName"] <- paste(MOM.Mat[,"Genus"],MOM.Mat[,"Species"], sep="_")
#MOM.Mat[,"TaxonomicName"] <- gsub(pattern = " ", replacement = "_", x = MOM.Mat[,"TaxonomicName"],)
MOM.Mat[,"CombinedMass(kg)"] <- MOM.Mat[,"Combined.Mass..g."]/1000
MOM.Mat[,"logMass.Kg"] <- log10(MOM.Mat[,"CombinedMass(kg)"])
MOM.Mat <- MOM.Mat[,c("Order","FAMILY","Genus","Species","TaxonomicName","CombinedMass(kg)","logMass.Kg")]
#combine same species if multiple entires(e.g. Muntiacus_reevesi has separate entries for mainland Asia and Insular)
MOM.Mat <- aggregate(MOM.Mat[,-c(1:5)], by = list(MOM.Mat$Order, MOM.Mat$FAMILY,MOM.Mat$Genus,
MOM.Mat$Species, MOM.Mat$TaxonomicName), mean)
colnames(MOM.Mat)[1:5] <- c("Order","FAMILY","Genus","Species","TaxonomicName")
ecoMat <- read.csv("~/Dropbox/CalcanealGearRatio/ExtantUngulate_BM_Eco_Mendoza_etal2007.csv")
ecoMat["TaxonomicName"] <- paste(ecoMat[,"Genus"],paste("_",ecoMat[,"Species"],sep=""),sep="")
#############################################################################################
###Analysis
#############################################################################################
getOneSpecimenAverage <- function(this.specimen, calcmeasure) {
this.measurements <- c("Length.Calc", "Length..medial.", "Length..lateral.", "Width..between.prox.condyles.", "Width..between.distal.condyles.navicular.facet.", "Length.to.Sustentacular.Process", "Length.to.Calcaneoastragular.facet", "Calcaneoastragular.facet.to.cuboid.facet")
colMeans(calcmeasure[calcmeasure$Catalog.Number== this.specimen, this.measurements], na.rm=TRUE)
}
getSpecimenAverages <- function(calcmeasure) {
specimen.vec <- sort(unique(calcmeasure$Catalog.Number))
n <- t(sapply(specimen.vec, function(x) calcmeasure[calcmeasure$Catalog.Number==x, c("Museum.Acronym", "Catalog.Number", "Genus", "Species")][1,]))
m <- data.frame(n, t(sapply(specimen.vec, getOneSpecimenAverage, calcmeasure)))
m
}
getSpeciesAverages <- function(measureMat) {
rownames(measureMat) <- NULL
species.vec <- sort(unique(measureMat$taxon))
this.measurements <- c("Length.Calc", "Length..medial.", "Length..lateral.", "Width..between.prox.condyles.", "Width..between.distal.condyles.navicular.facet.", "Length.to.Sustentacular.Process", "Length.to.Calcaneoastragular.facet", "Calcaneoastragular.facet.to.cuboid.facet")
tst <- species.vec[1]
this.species <- tst
spMat <-matrix(nrow=length(species.vec),ncol = length(this.measurements))
for(jj in seq(1, length(species.vec),1))
{
this.species <- species.vec[jj]
spMat[jj,] <- colMeans(measureMat[measureMat$taxon == this.species, this.measurements], na.rm=TRUE)
}
colnames(spMat) <- this.measurements
rownames(spMat) <- species.vec
spMat
}
m <- getSpecimenAverages(calcmeasure)
m$taxon <- paste(m$Genus, tolower(m$Species), sep="_")
#m$R2 <- m$Length.to.Calcaneoastragular.facet/m$Calcaneoastragular.facet.to.cuboid.facet #this is based on the figure caption
m$R2 <- m$Calcaneoastragular.facet.to.cuboid.facet/m$Length.to.Calcaneoastragular.facet
#m$R2 <- m$Length.to.Calcaneoastragular.facet/m$Length.Calc #this is based on the designation listed in the figure
m$R3 <- m$Length.Calc/m$Length.to.Sustentacular.Process
m$A1 <- rowMeans(m[,c("Length..medial.", "Length..lateral.")], na.rm=TRUE)/rowMeans(m[,c("Width..between.prox.condyles.", "Width..between.distal.condyles.navicular.facet.")], na.rm=TRUE)
m$habitat <- ecoMat$Habitat[match(m$taxon, ecoMat$TaxonomicName)]
m$BM <- MOM.Mat$logMass.Kg[match(m$taxon, MOM.Mat$TaxonomicName)]
#family
m$Family <- MOM.Mat$FAMILY[match(m$taxon, MOM.Mat$TaxonomicName)]
#order
#family
m$Order <- MOM.Mat$Order[match(m$taxon, MOM.Mat$TaxonomicName)]
##remove Perissodactyls (3 species)
m <- m[!m$Order == "Perissodactyla",]
##species level
m.sp <- as.data.frame(getSpeciesAverages(m))
m.sp$taxon <- row.names(m.sp)
#m$R2 <- m$Length.to.Calcaneoastragular.facet/m$Calcaneoastragular.facet.to.cuboid.facet #this is based on the figure caption
m.sp$R2 <- m.sp$Calcaneoastragular.facet.to.cuboid.facet/m.sp$Length.to.Calcaneoastragular.facet
#m$R2 <- m$Length.to.Calcaneoastragular.facet/m$Length.Calc #this is based on the designation listed in the figure
m.sp$R3 <- m.sp$Length.Calc/m.sp$Length.to.Sustentacular.Process
m.sp$A1 <- rowMeans(m.sp[,c("Length..medial.", "Length..lateral.")], na.rm=TRUE)/rowMeans(m.sp[,c("Width..between.prox.condyles.", "Width..between.distal.condyles.navicular.facet.")], na.rm=TRUE)
m.sp$habitat <- ecoMat$Habitat[match(m.sp$taxon, ecoMat$TaxonomicName)]
m.sp$BM <- MOM.Mat$logMass.Kg[match(m.sp$taxon, MOM.Mat$TaxonomicName)]
#reorder
m.sp <- m.sp[,c(9,13,1:8,14, 10:12)]
m.sp <- m.sp[complete.cases(m.sp),] #68 species
m <- m.sp
#############################################################################################habitat.colors <- c("dodgerblue", "chartreuse", "goldenrod1")Polly 2010 Figure 4c
# habitat.colors <- rainbow(length(unique(m$habitat)))
plot(m$R2, m$R3, xlim=c(0.35,0.65), ylim=c(1.1, 1.5), pch=21, cex=1, bg=habitat.colors[m$habitat], xlab = "Gear Ratio (R2)", ylab = "Gear Ratio (R3)")
# plot(m$R2, m$R3, pch=21, cex=1, bg=habitat.colors[m$habitat])
text(m$R2, m$R3, labels=m$taxon, pos=3, cex=0.3)
abline(lm(m$R3 ~ m$R2))
Polly 2010 Figure 6
plot(sort(m$R3), pch=21, bg=habitat.colors[m$habitat], type="n", ylim=c(1,1.5))
rect(xleft=-10, ybottom=mean(m$R3, na.rm=TRUE)-sd(m$R3, na.rm=TRUE), xright=150, ytop=mean(m$R3, na.rm=TRUE)+sd(m$R3, na.rm=TRUE), col=adjustcolor("black", alpha.f=0.25))
abline(h=mean(m$R3, na.rm=TRUE), lty=2)
points(sort(m$R3), pch=21, bg=habitat.colors[m$habitat[order(m$R3)]])
#text(seq_len(nrow(m)), sort(m$R3), labels=m$taxon[order(m$R3)], pos=3, cex=0.2)
#axis(1, at = seq_len(nrow(m)), labels = sort(m$R3)) make diagonal labels for each taxonPolly 2010 Figure 6 but for R2, instead of R3
plot(sort(m$R2), pch=21, bg=habitat.colors[m$habitat], type="n", ylab = "Gear Ratio (R2)")
rect(xleft=-10, ybottom=mean(m$R2, na.rm=TRUE)-sd(m$R2, na.rm=TRUE), xright=150, ytop=mean(m$R2, na.rm=TRUE)+sd(m$R2, na.rm=TRUE), col=adjustcolor("black", alpha.f=0.25))
abline(h=mean(m$R2, na.rm=TRUE), lty=2)
points(sort(m$R2), pch=21, bg=habitat.colors[m$habitat[order(m$R2)]])
text(seq_len(nrow(m)), sort(m$R2), labels=m$taxon[order(m$R2)], pos=3, cex=0.3)
Polly 2010 Figure 6 but for A1, instead of R3
plot(sort(m$A1), pch=21, bg=habitat.colors[m$habitat], type="n", ylab = "Area of Astragalus")
rect(xleft=-10, ybottom=mean(m$A1, na.rm=TRUE)-sd(m$A1, na.rm=TRUE), xright=150, ytop=mean(m$A1, na.rm=TRUE)+sd(m$A1, na.rm=TRUE), col=adjustcolor("black", alpha.f=0.25))
abline(h=mean(m$A1, na.rm=TRUE), lty=2)
points(sort(m$A1), pch=21, bg=habitat.colors[m$habitat[order(m$A1)]])
text(seq_len(nrow(m)), sort(m$A1), labels=m$taxon[order(m$A1)], pos=3, cex=0.3)
### Correlation of ratios with body mass
this.ratio <- m$R2
plot(m$BM, this.ratio, pch=21, bg=habitat.colors[m$habitat], ylab = "Gear Ratio (R2)", xlab = "Body mass log(kg)")
text(m$BM, this.ratio, labels=m$taxon, pos=3, cex=0.3)
abline(lm(this.ratio ~ m$BM))
cor.test(m$BM, this.ratio)##
## Pearson's product-moment correlation
##
## data: m$BM and this.ratio
## t = -0.47015, df = 66, p-value = 0.6398
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.2921755 0.1831732
## sample estimates:
## cor
## -0.0577751 this.ratio <- m$R3
plot(m$BM, this.ratio, pch=21, bg=habitat.colors[m$habitat], ylab = "Gear Ratio (R3)", xlab = "Body mass log(kg)")
text(m$BM, this.ratio, labels=m$taxon, pos=3, cex=0.3)
abline(lm(this.ratio ~ m$BM))
cor.test(m$BM, this.ratio)##
## Pearson's product-moment correlation
##
## data: m$BM and this.ratio
## t = -1.2026, df = 66, p-value = 0.2334
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.37187876 0.09531537
## sample estimates:
## cor
## -0.14643763-Habitat Hypothesis
PCA
library(ggfortify)## Loading required package: ggplot2dropCol <- c(1:3,11:14)
df <- log(m[,-dropCol])
habitat <- m$habitat
df <- scale(df)
df.pca <- prcomp(df)
summary(df.pca)## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6
## Standard deviation 2.6306 0.20277 0.13787 0.10559 0.07592 0.04685
## Proportion of Variance 0.9886 0.00587 0.00272 0.00159 0.00082 0.00031
## Cumulative Proportion 0.9886 0.99442 0.99714 0.99873 0.99955 0.99987
## PC7
## Standard deviation 0.03068
## Proportion of Variance 0.00013
## Cumulative Proportion 1.00000#plot(df.pca$x[,1], df.pca$x[,2], pch = 16, col = habitat.colors[m$habitat], ylab = "PC2", xlab = "PC1")
autoplot(df.pca, data = m, colour = 'habitat', label = FALSE, label.size = 3,
loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3)+
scale_color_manual(breaks = c("CH", "MH", "OH"),
values=c("dodgerblue", "chartreuse", "goldenrod1"))
#plot(df.pca$x[,2], df.pca$x[,3], pch = 16, col = habitat.colors[m$habitat], ylab = "PC3", xlab = "PC2")
autoplot(df.pca, data = m, x=2, y=3, choices = c(2:3), colour = 'habitat',
label = TRUE, label.size = 3,
loadings = TRUE, loadings.label = TRUE, loadings.label.size = 3) +
scale_color_manual(breaks = c("CH", "MH", "OH"),
values=c("dodgerblue", "chartreuse", "goldenrod1"))
K-means
set.seed(1)
#set-symbols for pch
ksym <- gsub("OH",16, m$habitat)
ksym <- gsub("MH",17, ksym)
ksym <- gsub("CH",18, ksym)
ksym <- as.numeric(ksym)
df.kmean <- kmeans(df, 3)
kcol <- gsub(1, "red", df.kmean$cluster)
kcol <- gsub(2, "darkgreen", kcol)
kcol<- gsub(3, "blue", kcol)
autoplot(df.kmean, data = m, size = 0.1) +
geom_point(shape = ksym, colour = kcol, size = 2) 
Discriminant Analysis
Linear Discriminant Analysis
library(lfda)
library(MASS)
lda.tst <- lda(x = as.matrix(df), grouping = m$habitat)
plot(lda.tst, col = habitat.colors[m$habitat])
Local Fisher Discriminant Analysis (LFDA)
model1_3h <- lfda(df, m$habitat, 3, metric="plain")
autoplot(model1_3h, data = m, frame = TRUE, frame.colour = 'habitat', size=0.1)+
geom_point(shape = ksym, size = 2) 
autoplot(model1_3h, data = m, x=2, y=3,frame = TRUE, frame.colour = 'habitat', size=0.1)+
geom_point(shape = ksym, size=2) 
###Kernal Local Fisher Discriminant Analysis (KLFDA)
model2_3h <- klfda(kmatrixGauss(df), m$habitat, 3, metric="plain")
autoplot(model2_3h, data = m, frame = TRUE, frame.colour = 'habitat', size=0.1)+
geom_point(shape = ksym, size = 2) 
autoplot(model2_3h, data = m,x=2,y=3, frame = TRUE, frame.colour = 'habitat', size=0.1)+
geom_point(shape = ksym, size = 2) 
Semi-supervised Local Fisher Discriminant Analysis (SELF)
model3_3h <- self(df, m$habitat, beta = 0.1, r = 3, metric="plain")
autoplot(model3_3h, data = m, frame = TRUE, frame.colour = 'habitat', size=0.1)+
geom_point(shape = ksym, size = 2) 
require(phytools)## Loading required package: phytools## Loading required package: ape## Loading required package: maps#read in tree
testTree <- read.nexus("/Users/emdoughty/Downloads/Appendix 3 alignment.nex")
testTreeTips <- testTree$tip.label
getTipTaxon <- function(treetips){
tipList <- strsplit(treetips,split = "_")
tipMat <- matrix(nrow = length(tipList), ncol = 2)
for(ii in seq(1, length(tipList),1))
{
tipMat[ii,1] <- tipList[[ii]][1]
tipMat[ii,2] <- tipList[[ii]][2]
}
labelTipsNew <- paste(tipMat[,1], tipMat[,2], sep="_")
return (labelTipsNew)
}
labelTipsNew <- getTipTaxon(testTree$tip.label)
#how many entries are present per species name
new.table <- table(labelTipsNew)
removeVec <- vector()
for(yy in seq(1, length(new.table),1))
{
if(new.table[yy] > 1)
{
names(new.table[yy])
tip.options <- which(labelTipsNew == names(new.table[yy]))
tip.remove <- sample(tip.options,size = length(tip.options)-1)
removeVec <- append(removeVec, tip.remove)
}
}
#check removeVec against names
testTree$tip.label[removeVec]## [1] "Alcelaphus_buselaphus_buselaphus"
## [2] "Balaenoptera_edeni_brydei_NC006928"
## [3] "Bos_javanicus_javanicus"
## [4] "Bos_taurus_taurus_EU177832"
## [5] "Bubalus_bubalis_bubalis_AF547270"
## [6] "Cephalophus_callipygus_1"
## [7] "Cephalophus_callipygus_2"
## [8] "Cervus_nippon_taiouanus_NC008462"
## [9] "Cervus_nippon_centralis_NC006993"
## [10] "Connochaetes_taurinus_2"
## [11] "Dama_dama_dama"
## [12] "Eudorcas_rufifrons_2"
## [13] "Gazella_dorcas_pelzelnii"
## [14] "Gazella_gazella_erlangeri"
## [15] "Gazella_subgutturosa_marica"
## [16] "Giraffa_camelopardalis_angolensis_NC012100"
## [17] "Mazama_americana_2"
## [18] "Mazama_nemorivaga_1"
## [19] "Naemorhedus_griseus_griseus_FJ207532"
## [20] "Odocoileus_virginianus_3"
## [21] "Odocoileus_virginianus_1"
## [22] "Pecari_tajacu_2_NC012103"
## [23] "Pecari_tajacu_3"
## [24] "Philantomba_monticola_2"
## [25] "Tragelaphus_scriptus_2"
## [26] "Tragelaphus_scriptus_3"#collapse/remove species with multiple entries
testTree_crop <- drop.tip(testTree, removeVec)
#check if it worked
newNameList <- getTipTaxon(testTree_crop$tip.label)
table(newNameList)## newNameList
## Addax_nasomaculatus Aepyceros_melampus
## 1 1
## Alcelaphus_buselaphus Alces_alces
## 1 1
## Ammotragus_lervia Antidorcas_marsupialis
## 1 1
## Antilocapra_americana Antilope_cervicapra
## 1 1
## Arabitragus_jayakari Axis_axis
## 1 1
## Axis_porcinus Balaena_mysticetus
## 1 1
## Balaenoptera_acutorostrata Balaenoptera_bonaerensis
## 1 1
## Balaenoptera_borealis Balaenoptera_edeni
## 1 1
## Balaenoptera_musculus Balaenoptera_omurai
## 1 1
## Balaenoptera_physalus Berardius_bairdii
## 1 1
## Bison_bison Bison_bonasus
## 1 1
## Blastocerus_dichotomus Bos_gaurus
## 1 1
## Bos_grunniens Bos_javanicus
## 1 1
## Bos_taurus Boselaphus_tragocamelus
## 1 1
## Bubalus_bubalis Bubalus_depressicornis
## 1 1
## Budorcas_taxicolor Camelus_bactrianus
## 1 1
## Camelus_dromedarius Camelus_ferus
## 1 1
## Caperea_marginata Capra_caucasica
## 1 1
## Capra_falconeri Capra_hircus
## 1 1
## Capra_ibex Capra_nubiana
## 1 1
## Capra_pyrenaica Capra_sibirica
## 1 1
## Capreolus_capreolus Capricornis_crispus
## 1 1
## Capricornis_milneedwardsii Capricornis_swinhoei
## 1 1
## Cephalophus_adersi Cephalophus_callipygus
## 1 1
## Cephalophus_dorsalis Cephalophus_jentinki
## 1 1
## Cephalophus_leucogaster Cephalophus_natalensis
## 1 1
## Cephalophus_nigrifrons Cephalophus_ogilbyi
## 1 1
## Cephalophus_rufilatus Cephalophus_silvicultor
## 1 1
## Cephalophus_spadix Cephalorhynchus_heavisidii
## 1 1
## Cervus_elaphus Cervus_nippon
## 1 1
## Choeropsis_liberiensis Connochaetes_gnou
## 1 1
## Connochaetes_taurinus Dama_dama
## 1 1
## Damaliscus_pygargus Delphinus_capensis
## 1 1
## Dorcatragus_megalotis Elaphodus_cephalophus
## 1 1
## Elaphurus_davidianus Equus_asinus
## 1 1
## Eschrichtius_robustus Eubalaena_australis
## 1 1
## Eubalaena_japonica Eudorcas_rufifrons
## 1 1
## Gazella_bennettii Gazella_cuvieri
## 1 1
## Gazella_dorcas Gazella_gazella
## 1 1
## Gazella_leptoceros Gazella_spekei
## 1 1
## Gazella_subgutturosa Giraffa_camelopardalis
## 1 1
## Grampus_griseus Hemitragus_jemlahicus
## 1 1
## Hippocamelus_antisensis Hippopotamus_amphibius
## 1 1
## Hippotragus_equinus Hippotragus_niger
## 1 1
## Hydropotes_inermis Hyemoschus_aquaticus
## 1 1
## Hyperoodon_ampullatus Inia_geoffrensis
## 1 1
## Kobus_ellipsiprymnus Kobus_leche
## 1 1
## Kogia_breviceps Lagenorhynchus_albirostris
## 1 1
## Lama_guanicoe Lipotes_vexillifer
## 1 1
## Litocranius_walleri Madoqua_kirkii
## 1 1
## Madoqua_saltiana Mazama_americana
## 1 1
## Mazama_gouazoubira Mazama_nemorivaga
## 1 1
## Mazama_rufina Megaptera_novaeangliae
## 1 1
## Monodon_monoceros Moschus_berezovskii
## 1 1
## Moschus_moschiferus Muntiacus_crinifrons
## 1 1
## Muntiacus_muntjak Muntiacus_reevesi
## 1 1
## Muntiacus_vuquangensis Naemorhedus_baileyi
## 1 1
## Naemorhedus_griseus Nanger_dama
## 1 1
## Nanger_granti Nanger_soemmerringii
## 1 1
## Neotragus_batesi Nesotragus_moschatus
## 1 1
## Odocoileus_hemionus Odocoileus_virginianus
## 1 1
## Okapia_johnstoni Oreamnos_americanus
## 1 1
## Oreotragus_oreotragus Oryx_beisa
## 1 1
## Oryx_dammah Oryx_gazella
## 1 1
## Oryx_leucoryx Ourebia_ourebi
## 1 1
## Ovibos_moschatus Ovis_aries
## 1 1
## Ozotoceros_bezoarcticus Panthera_tigris
## 1 1
## Pantholops_hodgsonii Pecari_tajacu
## 1 1
## Pelea_capreolus Phacochoerus_africanus
## 1 1
## Philantomba_maxwelli Philantomba_monticola
## 1 1
## Phocoena_phocoena Physeter_macrocephalus
## 1 1
## Platanista_gangetica Pontoporia_blainvillei
## 1 1
## Potamochoerus_porcus Procapra_gutturosa
## 1 1
## Przewalskium_albirostris Pseudois_nayaur
## 1 1
## Pseudoryx_nghetinhensis Pteropus_dasymallus
## 1 1
## Pudu_mephistophiles Pudu_puda
## 1 1
## Rahicerus_campestris Rangifer_tarandus
## 1 1
## Redunca_arundinum Redunca_fulvorufula
## 1 1
## Rucervus_duvauceli Rucervus_eldi
## 1 1
## Rupicapra_pyrenaica Rupicapra_rupicapra
## 1 1
## Rusa_alfredi Rusa_timorensis
## 1 1
## Rusa_unicolor Saiga_tatarica
## 1 1
## Sousa_chinensis Stenella_attenuata
## 1 1
## Stenella_coeruleoalba Sus_barbatus
## 1 1
## Sus_scrofa Sylvicapra_grimmia
## 1 1
## Syncerus_caffer Tetracerus_quadricornis
## 1 1
## Tragelaphus_angasii Tragelaphus_derbianus
## 1 1
## Tragelaphus_eurycerus Tragelaphus_imberbis
## 1 1
## Tragelaphus_oryx Tragelaphus_scriptus
## 1 1
## Tragelaphus_spekii Tragelaphus_strepsiceros
## 1 1
## Tragulus_kanchil Tursiops_aduncus
## 1 1
## Tursiops_truncatus Vicugna_pacos
## 1 1testTree_crop$tip.label <- newNameList
#find which species are not on the tree
unique(m$taxon[!m$taxon %in% newNameList])## [1] "Beatragus_hunteri" "Cephalophus_zebra"
## [3] "Cervus_canadensis" "Damaliscus_lunatus"
## [5] "Hemitragus_hylocrius" "Kobus_kob"
## [7] "Lama_pacos" "Madoqua_guentheri"
## [9] "Neotragus_moschatus" "Ovis_canadensis"
## [11] "Ovis_dalli" "Ovis_orientalis"
## [13] "Phacochoerus_aethiopicus" "Philantomba_maxwellii"
## [15] "Raphicerus_campestris" "Rucervus_eldii"
## [17] "Taurotragus_oryx" "Tayassu_pecari"
## [19] "Tragulus_javanicus" "Tragulus_napu"
## [21] "Vicugna_vicugna"unique(testTree_crop$tip.label[!testTree_crop$tip.label %in% m$taxon])## [1] "Pteropus_dasymallus" "Equus_asinus"
## [3] "Panthera_tigris" "Phacochoerus_africanus"
## [5] "Potamochoerus_porcus" "Sus_barbatus"
## [7] "Lama_guanicoe" "Vicugna_pacos"
## [9] "Camelus_dromedarius" "Camelus_ferus"
## [11] "Choeropsis_liberiensis" "Balaena_mysticetus"
## [13] "Eubalaena_australis" "Eubalaena_japonica"
## [15] "Caperea_marginata" "Balaenoptera_acutorostrata"
## [17] "Balaenoptera_bonaerensis" "Eschrichtius_robustus"
## [19] "Megaptera_novaeangliae" "Balaenoptera_physalus"
## [21] "Balaenoptera_musculus" "Balaenoptera_omurai"
## [23] "Balaenoptera_borealis" "Balaenoptera_edeni"
## [25] "Kogia_breviceps" "Physeter_macrocephalus"
## [27] "Platanista_gangetica" "Hyperoodon_ampullatus"
## [29] "Berardius_bairdii" "Lipotes_vexillifer"
## [31] "Inia_geoffrensis" "Pontoporia_blainvillei"
## [33] "Monodon_monoceros" "Phocoena_phocoena"
## [35] "Lagenorhynchus_albirostris" "Cephalorhynchus_heavisidii"
## [37] "Grampus_griseus" "Stenella_attenuata"
## [39] "Sousa_chinensis" "Tursiops_truncatus"
## [41] "Stenella_coeruleoalba" "Delphinus_capensis"
## [43] "Tursiops_aduncus" "Tragulus_kanchil"
## [45] "Okapia_johnstoni" "Capreolus_capreolus"
## [47] "Ozotoceros_bezoarcticus" "Mazama_gouazoubira"
## [49] "Hippocamelus_antisensis" "Blastocerus_dichotomus"
## [51] "Mazama_nemorivaga" "Mazama_rufina"
## [53] "Pudu_mephistophiles" "Elaphodus_cephalophus"
## [55] "Muntiacus_muntjak" "Muntiacus_crinifrons"
## [57] "Muntiacus_vuquangensis" "Rucervus_duvauceli"
## [59] "Elaphurus_davidianus" "Rucervus_eldi"
## [61] "Rusa_alfredi" "Cervus_elaphus"
## [63] "Przewalskium_albirostris" "Moschus_berezovskii"
## [65] "Moschus_moschiferus" "Tetracerus_quadricornis"
## [67] "Tragelaphus_scriptus" "Tragelaphus_angasii"
## [69] "Tragelaphus_oryx" "Tragelaphus_derbianus"
## [71] "Syncerus_caffer" "Bubalus_bubalis"
## [73] "Bubalus_depressicornis" "Pseudoryx_nghetinhensis"
## [75] "Bison_bonasus" "Bos_taurus"
## [77] "Bos_grunniens" "Bos_gaurus"
## [79] "Neotragus_batesi" "Nesotragus_moschatus"
## [81] "Procapra_gutturosa" "Rahicerus_campestris"
## [83] "Dorcatragus_megalotis" "Madoqua_saltiana"
## [85] "Madoqua_kirkii" "Litocranius_walleri"
## [87] "Eudorcas_rufifrons" "Nanger_granti"
## [89] "Nanger_dama" "Gazella_spekei"
## [91] "Gazella_dorcas" "Gazella_gazella"
## [93] "Gazella_bennettii" "Gazella_leptoceros"
## [95] "Gazella_cuvieri" "Redunca_fulvorufula"
## [97] "Redunca_arundinum" "Pelea_capreolus"
## [99] "Kobus_leche" "Oreotragus_oreotragus"
## [101] "Philantomba_maxwelli" "Sylvicapra_grimmia"
## [103] "Cephalophus_dorsalis" "Cephalophus_jentinki"
## [105] "Cephalophus_silvicultor" "Cephalophus_spadix"
## [107] "Cephalophus_adersi" "Cephalophus_ogilbyi"
## [109] "Cephalophus_callipygus" "Cephalophus_leucogaster"
## [111] "Cephalophus_natalensis" "Damaliscus_pygargus"
## [113] "Alcelaphus_buselaphus" "Connochaetes_gnou"
## [115] "Hippotragus_equinus" "Hippotragus_niger"
## [117] "Addax_nasomaculatus" "Oryx_beisa"
## [119] "Oryx_dammah" "Pantholops_hodgsonii"
## [121] "Capricornis_crispus" "Capricornis_milneedwardsii"
## [123] "Capricornis_swinhoei" "Naemorhedus_baileyi"
## [125] "Naemorhedus_griseus" "Ovis_aries"
## [127] "Ammotragus_lervia" "Arabitragus_jayakari"
## [129] "Rupicapra_pyrenaica" "Budorcas_taxicolor"
## [131] "Pseudois_nayaur" "Capra_sibirica"
## [133] "Capra_nubiana" "Capra_pyrenaica"
## [135] "Capra_ibex" "Capra_caucasica"
## [137] "Capra_falconeri"#remove from each dataset after removing incompletes in dataset
m <- m[complete.cases(m),]
m.tree <- m[m$taxon %in% newNameList,]
df.tree <- df[rownames(df) %in% newNameList,]
testTree_crop <- drop.tip(testTree_crop,tip = testTree_crop$tip.label[!testTree_crop$tip.label %in% m.tree$taxon])
artTree <- testTree_crop
#align so that tips and m.tree are in same order
m.tree <- m.tree[artTree$tip.label,]
df.tree <- df.tree[artTree$tip.label,]plot(artTree, cex = 0.4)
df.phylo.pca <- phyl.pca(artTree, df.tree, method="BM", mode="cov")
summary(df.phylo.pca)## Importance of components:
## PC1 PC2 PC3 PC4
## Standard deviation 6.0083652 0.63086636 0.389308403 0.251244557
## Proportion of Variance 0.9816734 0.01082254 0.004121373 0.001716515
## Cumulative Proportion 0.9816734 0.99249599 0.996617358 0.998333874
## PC5 PC6 PC7
## Standard deviation 0.1896223710 0.1385120336 0.078285155
## Proportion of Variance 0.0009777629 0.0005217103 0.000166653
## Cumulative Proportion 0.9993116367 0.9998333470 1.000000000plot(df.phylo.pca)
#plot(df.phylo.pca$S[,1],df.phylo.pca$S[,2], col = habitat.colors[m.tree$habitat],
# pch =16)
#text(df.phylo.pca$S[,1],df.phylo.pca$S[,2], labels=artTree$tip.label, cex = 0.5)require(maptools)## Loading required package: maptools## Loading required package: sp## Checking rgeos availability: FALSE
## Note: when rgeos is not available, polygon geometry computations in maptools depend on gpclib,
## which has a restricted licence. It is disabled by default;
## to enable gpclib, type gpclibPermit()plot(df.phylo.pca$S[,1],df.phylo.pca$S[,2], col = habitat.colors[m.tree$habitat],
pch =16, xlab = "PC1 (98.73%)", ylab = "PC2 (0.67 %)")
#points(df.phylo.pca$L[,1:2])
#pointLabel(df.phylo.pca$L[,2:3], # set position of labels
# labels=rownames(df.phylo.pca$L), # print labels
# cex=0.75 )plot(df.phylo.pca$S[,2],df.phylo.pca$S[,3], col = habitat.colors[m.tree$habitat],
pch =16, xlab = "PC2 (0.67%)", ylab = "PC3 (0.45%)")
#points(df.phylo.pca$L[,2:3])
#pointLabel(df.phylo.pca$L[,2:3], # set position of labels
# labels=rownames(df.phylo.pca$L), # print labels
# cex=0.75 )
#pointLabel(df.phylo.pca$S[,2:3], # set position of labels
# labels=rownames(df.phylo.pca$S), # print labels
# cex=0.75)#loadings
#quartz(height = 6, width = 10)
par(mfrow=c(2,2))
plot(df.phylo.pca$L[,1],df.phylo.pca$L[,2], xlab="PC1", ylab="PC2", main = "Phylogenetic PCA")
pointLabel(df.phylo.pca$L[,1:2], # set position of labels
labels=rownames(df.phylo.pca$L), # print labels
cex=0.75 )
plot(df.phylo.pca$L[,2],df.phylo.pca$L[,3],xlab="PC2", ylab="PC3")
pointLabel(df.phylo.pca$L[,2:3], # set position of labels
labels=rownames(df.phylo.pca$L), # print labels
cex=0.75 )